/*
 * Copyright 2010-2012 Susanta Tewari. <freecode4susant@users.sourceforge.net>
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

package provider.impl.genomemap.data;

import commons.core.Simulator;
import commons.provider.ProviderException;
import commons.provider.impl.BaseProviderImpl;
import genomemap.data.TAB08Data;
import genomemap.model.TAB08;
import genomemap.provider.Providers;
import org.apache.commons.lang.ArrayUtils;

import java.util.ArrayList;
import java.util.List;

/**
 * We simulate data using a different poisson process in each interval.
 * We want to check if we can estimate the probability of a crossover(a pseudo-crossover)
 * closely to the model parameter.Note that occurrence of a single crossover is
 * not the only way a recombinant can be produced, it depends on the strands involved
 * in the crossover. For example, if all the four strands are involved in 2 crossovers
 * then all four meiotic products are recombinant. In other words, occurrence
 * of 2 crossovers does not guarantee that all products will be parental. The interpretation
 * of the model parameter 'c' is  thus the probability of a chiasmata in the corresponding
 * interval. Note that significantly different values of c would indicate that
 * chiasmata probabilities change across intervals.
 *
 * @version 1.3.0
 * @author Susanta Tewari
 * @history Created on 1/7/13
 * @since 1.3.0
 */
public class TAB08DataSimulatorProviderImpl_1
        extends BaseProviderImpl<Simulator.DATA<TAB08, TAB08Data>>
        implements Providers.TAB08DataSimulatorProvider {

    /** sample size */
    private int n;

    /** exchange probability in each interval */
    private double[] cProb;

    public TAB08DataSimulatorProviderImpl_1() {}

    /**
     *
     * @param input
     * @param sampleSize
     */
    TAB08DataSimulatorProviderImpl_1(final TAB08 input, final int sampleSize) {

        n     = sampleSize;
        cProb = ArrayUtils.toPrimitive(input.getExchangeProbs().toArray(new Double[0]));
    }

    /**
     * data: containing genotypes. Rows indicate meioses and columns are for genes
     * tetrad : 4 by loci created in a meiosis
     * @return simulated data
     */
    private byte[][] simulate() {

        final int loci      = cProb.length + 1;
        final byte[][] data = new byte[n][loci];
        byte[][] tetrad;
        int[] iAndj;

        for (int i = 0; i < n; i++) {

            tetrad = new byte[4][loci];


            // initialize tetrad the first column
            tetrad[0][0] = 1;
            tetrad[1][0] = 1;
            tetrad[2][0] = 0;
            tetrad[3][0] = 0;

            for (int index = 1; index < loci; index++) {

                iAndj = getIandJ(cProb[index - 1]);

                switch (iAndj[0]) {

                    case 0 :
                        tetrad[0][index] = tetrad[0][index - 1];
                        tetrad[1][index] = tetrad[1][index - 1];
                        tetrad[2][index] = tetrad[2][index - 1];
                        tetrad[3][index] = tetrad[3][index - 1];

                        break;

                    case 1 :
                        tetrad[0][index] = tetrad[2][index - 1];
                        tetrad[1][index] = tetrad[1][index - 1];
                        tetrad[2][index] = tetrad[0][index - 1];
                        tetrad[3][index] = tetrad[3][index - 1];

                        break;

                    case 2 :
                        tetrad[0][index] = tetrad[0][index - 1];
                        tetrad[1][index] = tetrad[2][index - 1];
                        tetrad[2][index] = tetrad[1][index - 1];
                        tetrad[3][index] = tetrad[3][index - 1];

                        break;

                    case 3 :
                        tetrad[0][index] = tetrad[0][index - 1];
                        tetrad[1][index] = tetrad[3][index - 1];
                        tetrad[2][index] = tetrad[2][index - 1];
                        tetrad[3][index] = tetrad[1][index - 1];

                        break;

                    case 4 :
                        tetrad[0][index] = tetrad[3][index - 1];
                        tetrad[1][index] = tetrad[1][index - 1];
                        tetrad[2][index] = tetrad[2][index - 1];
                        tetrad[3][index] = tetrad[0][index - 1];

                        break;
                }

                byte swap = 0;

                switch (iAndj[1]) {

                    case 0 :
                        tetrad[0][index] = tetrad[0][index];
                        tetrad[1][index] = tetrad[1][index];
                        tetrad[2][index] = tetrad[2][index];
                        tetrad[3][index] = tetrad[3][index];

                        break;

                    case 1 :
                        swap             = tetrad[0][index];
                        tetrad[0][index] = tetrad[2][index];
                        tetrad[1][index] = tetrad[1][index];
                        tetrad[2][index] = swap;
                        tetrad[3][index] = tetrad[3][index];

                        break;

                    case 2 :
                        swap             = tetrad[1][index];
                        tetrad[0][index] = tetrad[0][index];
                        tetrad[1][index] = tetrad[2][index];
                        tetrad[2][index] = swap;
                        tetrad[3][index] = tetrad[3][index];

                        break;

                    case 3 :
                        swap             = tetrad[1][index];
                        tetrad[0][index] = tetrad[0][index];
                        tetrad[1][index] = tetrad[3][index];
                        tetrad[2][index] = tetrad[2][index];
                        tetrad[3][index] = swap;

                        break;

                    case 4 :
                        swap             = tetrad[0][index];
                        tetrad[0][index] = tetrad[3][index];
                        tetrad[1][index] = tetrad[1][index];
                        tetrad[2][index] = tetrad[2][index];
                        tetrad[3][index] = swap;

                        break;
                }
            }

            data[i] = tetrad[randomIndex()];
        }

        return data;
    }

    /**
     * Method description
     *
     * @param c description
     *
     * @return description
     */
    private static int[] getIandJ(final double c) {


        // iAndj = i , j could be any one of values from 0 to 4.
        final int[] result = new int[2];
        double rnd         = Math.random();

        if (rnd <= (1 - c)) {
            result[0] = 0;
        } else {


            // chosse a strand out of 4 randomly for first crossover
            rnd = Math.random();

            if (rnd <= 0.25) result[0] = 1;
            if ((rnd <= 0.50) && (rnd > 0.25)) result[0] = 2;
            if ((rnd <= 0.75) && (rnd > 0.50)) result[0] = 3;
            if (rnd > 0.75) result[0] = 4;
        }

        rnd = Math.random();

        if (rnd <= (1 - c)) {
            result[1] = 0;
        } else {


            // chosse a strand out of 4 randomly for first crossover
            rnd = Math.random();

            if (rnd <= 0.25) result[1] = 1;
            if ((rnd <= 0.50) && (rnd > 0.25)) result[1] = 2;
            if ((rnd <= 0.75) && (rnd > 0.50)) result[1] = 3;
            if (rnd > 0.75) result[1] = 4;
        }

        return result;
    }

    /**
     * randomIndex(): gives one of the values of 0,1,2 and 3 equally likely,using
     * no chromatid inteference.
     *
     * @return description
     */
    private static int randomIndex() {

        int result       = 0;
        final double rnd = Math.random();

        if (rnd <= 0.25) result = 0;
        if ((rnd <= 0.50) && (rnd > 0.25)) result = 1;
        if ((rnd <= 0.75) && (rnd > 0.50)) result = 2;
        if (rnd > 0.75) result = 3;

        return result;
    }

    @Override
    public Simulator.DATA<TAB08, TAB08Data> create() throws ProviderException {

        return new Simulator.DATA<TAB08, TAB08Data>() {

            @Override
            public Iterable<TAB08Data> simulate(final DataInput<TAB08> input, final int simCount) {

                final List<TAB08Data> result = new ArrayList<>();
                final TAB08 model            = input.getModel();
                final TAB08DataSimulatorProviderImpl_1 simulator =
                    new TAB08DataSimulatorProviderImpl_1(model, input.getSampleSize());

                for (int i = 0; i < simCount; i++) {
                    result.add(new TAB08Data(model, simulator.simulate()));
                }

                return result;
            }
        };
    }

    @Override
    public Class<TAB08Data> getDataClass() {
        return null;
    }

    @Override
    public Class<TAB08> getModelClass() {
        return TAB08.class;
    }
}
